CHAPITRE 1 


ANALYSE D’ERREURS 


1.1 REPRESENTATION DES NOMBRES DANS UN CALCULATEUR 

Un calculateur ne peut fournir que des réponses approximatives. 

Les approximations utilisées dépendent à la fois des contraintes physiques (espace mémoire 
...) et du choix des méthodes retenues par le créateur du programme. 

La première contrainte est que le système numérique d’un calculateur quelconque est 
discret. C’est-à-dire qu’il ne comporte qu’un nombre fini de nombres. 

Il en découle que, sauf dans les cas les plus simples, tous les calculs seront entachés 
d’erreurs. 

Exemple 

Nous voulons calculer la valeur de (llllllll) 2 . 

Sur une calculatrice, on trouve une valeur approchée comme 1.234567876543e + 14 
La valeur exacte est 123456787654321 

1.1.1 Stockage des nombres 

Les nombres sont stockés dans un ordinateur comme entiers ou réels. 

Les nombres entiers 

Les nombres entiers sont ceux auquels nous sommes habitués sauf que le plus grand 
nombre représentable dépend des nombres d’octets utilisés. 

- Avec deux octets, on peut représenter les entiers compris entre —32768 et 32767 

- Avec quatre octets, on peut représenter les entiers compris entre —2147483648 et 
2147483647 
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Les nombres réels 

Dans la mémoire d’un calculateur les nombres réels sont représentés en notation 
flottante. 

En notation flottante, un nombre a la forme : 

a — Y x ¥ 

. b est la base du système numérique utilisé. 

. Y est la mantisse. 

. e est l’exposant. 

La base 

La base est celle dans laquelle on écrit les nombres. (La plupart des calculatrices 
utilisent le système décimal c’est-à-dire b = 10. Les ordinateurs utilisent le système 
binaire c’est-à-dire b = 2.) 

La mantisse 

La mantisse Y est un nombre de la forme : 

Y = ±0.did2 ■ ■ ■ dt 

où 0 < di < b et t est fixé et fini avec d\ ^ 0. 

Les di sont les chiffres décimaux de Y (ou les digits). 

11 faut, en général, un nombre infini de digits pour écrire exactement un nombre 
réel quelconque (on en a calculé des millions pour fl). 

Dans une calculatrice, le nombre t est généralement voisin de 10. 

Dans les grands ordinateurs, ce nombre peut prendre habituellement deux valeurs, 
la plus petite correspond à ce que l’on appelle la simple précision et la plus grande 
à la double précision. 

- L’exposant 

L’exposant e donne l’ordre de grandeur du nombre. On a : 

m < e < M 

où m et M sont des caractéristiques du calculateur. 
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Si e > M on dit qu’il y a dépassement (ou surpassement) de capacité. 

Si e < m on dit qu’il y a sous passement de capacité. 

Dans le cas d’un dépassement de capacité dans l’ordinateur il y a arrêt des calculs 
et l’impression d’un message d’erreur. 

Dans le cas d’un sous passement de capacité certains ordinateurs s’arrêtent et 
imprime un message d’erreur tandis que d’autres remplacent le nombre en cause 
par zéro et continuent les calculs. 

Si x est un nombre fourni à un ordinateur. On note fl(x ) sa représentation en virgule 
flottante. 

Troncature d’un nombre 

Considérons le nombre x = A = 0.06666... 

15 

Nous aurons si, t — 5 

fl(x ) = 0.66666 x 10" 1 

Nous avons tous simplement négligé les décimaux que nous ne pouvions stocker. On 
dit que l’on tronque et l’opération s’appelle troncature. 

Par contre, dans les calculatrices, il y a la plupart du temps une opération d’arron- 
dissement. 

Arrondissement d’un nombre 

Dans une représentation arrondie, lorsque la première décimale négligée est supérieure 
ou égale à 5, on ajoute 1 à la dernière décimale conservée. 

Ainsi, pour x = A = 0.06666... nous aurons, si t = 5 

fl(x ) = 0.66667 x 10" 1 
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1.1.2 Erreur d’affectation 

Théorème 

Dans une arithmétique flottante à t digits on a : 

| x — fl(x) | < 5 | x | p lCT f 

avec p = 1 dans le cas de l’arrondi et p = 2 dans le cas de la troncature. 

1.2 

OPERATIONS ARITHMETIQUES ELEMENTAIRES EN VIRGULE FLOTTANTE 


L’addition et la soustraction flottante 

- Xi ® x 2 = fl(fl(xi) + fl(x 2 )) 

- Xi © x 2 = fl(fl(xi) - fl(x 2 )) 

Considérons X\ et x 2 tels que : | X\ | > | x 2 \ 
on a : 

fl( x i) = 10 ei Ui 

fl(x 2 ) = 10 e2 Y 2 = 10 ei Y' avec E 2 ' = 10 e2_ei Y 2 

Si les exposants ne sont pas les mêmes, on doit aligner, c’est-à-dire rendre le plus 
petit exposant égal au plus grand. 

Exemples 

On considère le cas t = 4 

1. fl(x i) = 10 5 (0.4316) fl(x 2 ) = 10" 1 (0.1852) 

On écrit : fl(x 2 ) = 10 5 (0.0000001852) 

fl(x i) + fl(x 2 ) = 10 5 (0.4316001852) 

D’où : X\ © x 2 = 10 5 (0.4316) en arrondissant ou en tronquant 

2. fl(x i) = 10 5 (0.4316) fl(x 2 ) = 10 2 (0.3422) 
fl(x i) + fl(x 2 ) = 10 5 (0.4319422) 

D’où : X\ © x 2 = 10 5 (0.4319) en arrondissant ou en tronquant 
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3. fl(x i) = 10 5 (0.4316) fl(x 2 ) = 10 5 (0.7511) 
fl(x i) + fl(x 2 ) = 10 5 (1.1827) = 10 6 (0.11827) 

D’où : Xi © x 2 = 10 6 (0.1182) par troncature 

= 10 6 (0.1183) par arrondi 

4. fl(x i) = 10 5 (0.4316) fl(x 2 ) = 10 5 (0.4315) 

fl(xi) ~ fl(x 2 ) = 10 5 (0.0001) 

D’où : x\ © x 2 = 10 2 (0.1000) en arrondissant ou en tronquant 

La multiplication flottante 

Xi ® X 2 = fl(fl(Xi) X fl(x 2 )) 

On a : 

fl(xi) x fl(x 2 ) = 10 ei+e2 l'i y 2 

Remarquons que 0.1 < | Di | , | Y 2 \ < 1 et que 0.01 < | Y\ Y 2 \ < 1 

Il sera par conséquent nécessaire, dans certains cas, de renormaliser la mantisse Y\ Y 2 
afin que son premier digit soit non nul. 

Exemples 

1. fl(x i) = 10 2 (0.2432) fl(x 2 ) = 10 3 (0.2000) 

fl(x i) x fl(x 2 ) = 10 5 (0.04864000) 

D’où : X\ ® x 2 — 10 4 (0.4864) par troncature ou arrondi 

2. fl(x i) = 10 2 (0.2432) fl(x 2 ) = 10 3 (0.6808) 

fl(x i) x fl(x 2 ) = 10 5 (0.16557056) 

D’où : x\ ® x 2 — 10 5 (0.1655) par troncature 

= 10 5 (0.1656) par arrondi 
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La division flottante 

Xi 0 x 2 = fl(fl(xi)/fl(x 2 )) 


On a : 

fK x l) _ \Ç\ e l~ e 2 
fK x 2 ) 


Ll 

x 2 


- si | Y 1 | < | Y 2 | alors 0.1 < ] || ] < 1. On calcule donc le quotient que l’on 
tronque ou l’on arrondit à t digits. L’exposant e est égal à e\ — e 2 - 

- si | | > | Y 2 | alors 1 < | j < 10. On remplace alors Y\ par Y[ = 10 —1 L 1 , on 

yt 

calcule yr et l’on tronque ou l’on arrondi ensuite à t digits. L’exposant e est égal à 

Cl — Ê2 + 1. 


Exemples 

1. fl(x 1 ) = 10 6 (0.4323) fl(x 2 ) = 10 5 (0.2000) 

= 10 1 (2.1615) = 10 2 (0.21615) 

D’où : x\ 0 x 2 = 10 2 (0.2161) par troncature 

= 10 2 (0.2162) par arrondi 

2. fl(x 1 ) = 10 6 (0.2539) fl(x 2 ) = 10 7 (0.3000) 

|g = 10- 1 (0.8463333...) 

D’où : x\ 0 x 2 = 10 _1 (0.8463) par troncature ou par arrondi 


Associativité de l’addition 

x + (y + z) peut être différent de (x + y) + z Soit, par exemple, à calculer la somme : 


1 + 0.0004 + 0.0006 = 1.001 

Avec un ordinateur pour lequel t = 4 en procédant par troncature. On a : 

1 © 0.0004 = 1 
(1 © 0.0004) © 0.0006 = 1 

0.0004 ©0.0006 = 0.001 
1© (0.0004 ©0.0006) = 1.001 

Cet exemple montre que l’addition flottante peut influencer le résultat de la sommation 
des séries à termes positifs. 


6 



Brahim Benouahmane 


Calcul de sommes à termes positifs, 
ordre de sommation 

On veut calculer : 


n 

S = a* a t > 0 

1=1 

En arithmétique finie, on calcule cette somme en formant la suite des sommes par- 
tielles : 


Si = //(ai) 
S 2 = fil © 02 


Sfc Sk — 1 O afc k 2, • • • , n 

Le résultat est alors ,5/ ~ S. 

L’ordre dans lequel on somme les a* peut changer la valeur de la somme S n car l’arith- 
métique flottante n’est pas associative. 

Exemple 


Soit 


n 


2=1 


+ Z 


n + 1 


Calculer cette somme, en arithmétique flottante, de deux façons différentes : 


Sn,l 

= i + - + -- 

• + 9 

2 

rr + n 

Sn, 2 

1 

1 „ 

— 9 + • 

• • + — + 1 

n 2 + n 

2 


avec n — 9 , 99 , 999 , 9999 


n 

Sn,l 

Sji,2 

Valeur exacte S 

9 

1.9000000000 

1.9000000000 

1.9 

99 

1.9900000000 

1.9900000000 

1.99 

999 

1.9990000000 

1.9990000000 

1.999 

9999 

1.9998999999 

1.9999000000 

1.9999 
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On voit que l’on obtient des résultats différents selon que l’on somme de 1 à n ou de 
n à 1 ;les meilleurs résultats étant obtenus dans le second cas. 

Perte de chiffres significatifs dans la soustraction 

Xi © X 2 = fl(fl(xi) - fl(x 2 )) 

Exemple 

Considérons les nombres ^7001 et \/7000. 

En arithmétique flottante à 8 chiffres, on a : 

V7ÔÔ1 ~ 0.83671979 x 10 2 
v 7 ™) ~ 0.83666003 x 10 2 

Donc 

V7ÔÔÏ - v 7 ™) = //((0. 83671979 - 0.83666003) x 10 2 ) = 0.59760000 x 10~ 2 
On peut obtenir un résultat plus précis en utilisant l’identité suivante : 

\fx ~ y/ÿ = (Vx -y/ÿ)x ^ 

On obtient alors : 

1 0 (\/7C)C)ï © V7ÔÜÜ) = 1 0 (0.16733798 x 10 3 ) 

= 0.59759297 x 10" 2 

La valeur exacte est 0.597592962824 x 10~ 2 

La soustraction est l’opération la plus dangereuse en calcul numérique. Elle peut am- 
plifier l’erreur relative de façon catastrophique. 

1.3 INSTABILITE NUMERIQUE 


Définition 
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On dira que le calcul ou l’algorithme est numériquement instable si de petits chan- 
gements dans les données entraînent de petits changements dans les résultats. Evidem- 
ment, dans le cas contraire on dira qu’il y a instabilité numérique. 

Exemple d’instabilité numérique 

On veut calculer 


f(n) = 


1,1 x n 


dx 


J o a + x 

Nous allons exprimer f(n) récursivement : 


a = cte > 1 


= C (f + « - a) ix = f x u-, dx _ a f ETL dx 

J o a + x J 0 J o a + x 


= af(n - 1) 

n 


n > 1 


/(O) = Ln(i^) 

L’algorithme fourni par cette relation est numériquement instable. 

Voici les résultats obtenus pour a = 10 et 
n — 0, 1, • • • ,12 


n 

f(n) calculé 

f(n) exact 

0 

0.0953102 

0.0953102 

1 

0.0468982 

0.0468982 

2 

0.0310180 

0.0310180 

3 

0.0231535 

0.0231535 

4 

0.0184647 

0.0184647 

5 

0.0153527 

0.0153529 

6 

0.0131401 

0.0131377 

7 

0.0114558 

0.0114806 

8 

0.0104421 

0.0101944 

9 

0.0066903 

0.0091672 

10 

0.0330968 

0.00832797 

11 

0.2400592 

0.00762944 

12 

2.4839249 

0.00703898 


à partir de n = 5, les valeurs calculées sont de moins en moins précises à chaque 
itération ; pour n > 10 les résultats obtenus sont complètement erronés. 
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Cet algorithme est d’autant plus instable que a est plus grand que 1. 

Pour ce faire supposons que l’erreur d’arrondi sur /(O) est égale à £ 0 et qu’aucune erreur 
n’est introduite dans les calculs subséquats. 

Notons /(n) les valeurs calculées. 

/(O) = /(O) + £ o 

f(n) = £ - a f(n - 1) n = 1,2, - -- 

Par suite, si r n désigne l’erreur sur f{n) 

r n = fin) - fin) = -a f(n - 1) + - - - + a f(n - 1) 

n n 

= -a(f(n- 1) -f{n- 1)) 

= -ar n _ i n = 1,2,--- 

et donc, puisque r 0 = nous trouvons r n = (— a) n Eq 

n = 1,2,--- 

L’erreur initiale est multipliée par un facteur a à chaque itération. 
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CHAPITRE 2 


INTERPOLATION POLYNOMIALE 


La façon la plus simple d’approcher une fonction est l’utilisation d’un polynôme d’in- 
terpolation. 

Le procédé est le suivant : 

1. On choisit n + 1 points distincts x 0 , Xi, • • • , x n 

2. On calcule y 0 = f(x 0 ),y 1 = f(xi), • • • ,y n = f(x n ) 

3. On cherche un polynôme P n de degré n tel que : P n (xi) = Vi , i — 0, 1, • • • ,n 

Nous allons montrer l’existence et l’unicité d’un tel polynôme en le construisant effec- 
tivement. 


2.1 Interpolation de Lagrange 


Théorème et définition 

Il existe un polynôme P n de degré n et un seul, tel que : 

P n (xi) = f(xi ) Vf = 0, 1, • • • , n 

n 

Ce polynôme s’écrit : P n (x) = E f{xi) Li(x ) 

2=0 

n 

avec LAx) = I [ — (polynômes de Lagrange) 

A. A. rp . rp . 

n 

Le polynôme P n s’appelle le polynôme d’interpolation de Lagrange de la fonc- 
tion / relativement aux points x 0 , Xi, ■ ■ ■ , x n . 

Preuve. 
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unicité : 

Supposons qu’il existe un polynôme P* de degré n vérifiant : P*(xj) = f(xi ) Vi = 
0, 1, - - - ,n 

Posons : Q n = P n - P*. 

Q n serait un polynôme de degré n vérifiant : Q n {xj) = 0 Vi = 0, 1, • • • , n de sorte 
que Q n aurait au moins (n+ 1) racines distinctes. 

Ceci n’est possible que si Q = 0 (polynôme identiquement nul) ce qui prouve l’unicité 
de P n . 

- existence : 

Nous voyons immédiatement que : L t est un polynôme de degré n et Li(xj ) = Sij 
(symbole de Kronecker). 

Il en résulte que le polynôme P n est de degré n et vérifie : P n (xi ) = f(xi) Vi = 
0, 1, - - - ,n 

Remarque 

Les (■ n + 1) polynômes de Lagrange sont linéairement indépendants et forment donc 
une base de l’espace vectoriel des polynômes de degré inférieur ou égal à n, appelée base 
de Lagrange. 

Calcul de P n (ot ) 

Le calcul de P n (a ) nécessite celui des (■ n + 1) quantités L^a), ce qui est coûteux. 

La méthode de Lagrange est d’un intérêt plus théorique que pratique car : 

- coûteux en nombres d’opérations 

- formulation peu aisée : si l’on ajoute un point x n+ \ les L, t doivent entièrement être 
recalculés. 


2.2 Erreur d’interpolation 

Etude de l’erreur e n (x) = f(x) — P n (x) (appelée erreur d’interpolation) pour tout 
x G [a, b] 
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Théorème 

Si /G C^ n+1 \[a, b]) alors pour tout x G [a, b], il existe £ x appartenant au plus petit 
intervalle fermé / contenant x, x 0 , Xi, ■ ■ ■ ,x n tel que : 

e n (x) = f{x) - P n (x) = F(x) (*) 

où F(x) = (x — x 0 ) (x — xi) ■ ■ ■ (x — x n ) 


Preuve. 

Si x = Xi alors e n (x) = 0 et l’égalité (*) est vérifiée trivialement ; 

Supposons que x ^ aq Vf = 0, 1, • • • , n et considérons pour x fixé la fonction g définie 
par : 

g(t) = e n (t) - Y^e n (x) 

La fonction g G C^ n+1 ^([a, b]) et s’annule en (n + 2) points distincts x, xo, xi, • • • , x n . 
Le théorème de Rolle montre que g' admet au moins (n + 1) racines dans I. 

D’où, en procédant par récurrence sur l’ordre de dérivation de g , la fonction r/ n+ L admet 
au moins une racine dans I. Soit £ x cette racine. On a : 

0 = 5 ,(n+1) 0vr) = / (n+1) (£x) - [j ^e n {x) 

D’où : e n (x) = F(x) 

2.3 Interpolation d’Hermite 

Problème : Déterminer un polynôme qui coïncide avec /, ainsi que sa dérivée avec /', 
aux points x 0 , Xi, ■ ■ ■ , x n . 

Théorème et définition 

Etant donnée une fonction / définie sur [a, b] et admettant des dérivées aux points x tl 
il existe un polynôme P -2 n + 1 de degré 2n + 1 et un seul tel que 

P 2 n+l(Xi) = f(Xi ) et P 2 n +1 (Xi) = ffa ) 
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Vi = 0, 1, • • • , n 

Le polynôme P2 n +i ainsi défini est appelé polynôme d’interpolation d’Hermite 

de la fonction / relativement aux points x 0 ,Xi, ■ ■ ■ ,x n . 

2.4 Erreur d’interpolation 

Etude de l’erreur e n (x) = f(x) — P n (x) pour tout x G [a, b] 

Théorème 

On considère le polynôme d’interpolation d’Hermite P-2n+\- 
On suppose que / G C (2n+2 ^([a, 6]) alors pour tout x G [a, b], il existe appartenant au 
plus petit intervalle fermé / contenant x, Xq, X\, ■ ■ ■ , x n tel que 

e n (x) =f(x)~ P 2n + iO) = F 2 (x) /( ( 2 ’^ ) 2 ( ) t ) 

Preuve. 

On considère pour tout x fixé distinct des Xi la fonction g de la variable t défini par 
g{t) = e n (t) - |^e n (x) 

La fonction g' G C (2n+1 ^([a, b]) et admet au moins (2 n + 2) zéros distincts. 

Le théorème de Rolle montre que g" admet au moins (2 n+ 1) racines dans I. 

D’où, en procédant par récurrence sur l’ordre de dérivation de g, la fonction g( 2n + 2 ) admet 
au moins une racine dans /. Soit £ x cette racine. On a : 

0 = J< 2 " +2| (&) = / <2 “ +2) (6) - ^ e n (x) 

D’où : e n (x) = F 2 (x) /< pL > 2 ( )t ) 

2.5 Interpolation itérée 

La détermination et l’évaluation du polynôme de Lagrange sont assez coûteuses lorsque 
le nombre de noeuds s’accroît. 

Nous allons développer un procédé itératif permettant de calculer le polynôme d’interpo- 
lation P n (x) basé sur n + 1 noeuds Xo, xi, ■ ■ ■ , x n à partir du polynôme P n -\(x) basé sur 
n noeuds xq, xi, ■ ■ ■ , x n -\. 
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Pour n > 1, le polynôme P n (x ) ~Pn-l (. x ) est de degré n et s’annule aux points 
x 0 , Xi, ■ ■ ■ , x n _i, il est donc de la forme : 

P n (x) - P n -i(x) = a n (x - x 0 ) (x - xi) • • • (x - x n _i) 

où a n est le coefficient dominant de P n {x) 

D’après la formule de Lagrange on a : 

n 

Pn(x) = y^J(x k ) L k (x) 

k = 0 

L k (x ) est un polynôme de degré n dont le coefficient dominant est : 

1 

{Xk-X o) -(x k -x k -l) (x k -x k+1 ) — (x k -x„) 


donc 


Qj n 


n 


Y. 


f(Xk) 

(x k -x 0 ) ■■■ (x k - x fc _i) (x k - x k+1 ) ■■■(x k - x n ) 


a 


n 


Le nombre a n est appelé différence divisée d’ordre n de la fonction f(x) et est noté : 
= f[x 0 ,xi, • • • ,x n ; 


Définissant f(x o) = f[x o], nous avons : 

P n (x) = P n -i(x) + (x-x 0 ) • • • ( x - x n _i) f[x o, x u *f , x n \ 

En utilisant ces relations pour n — 1, 2, • • • , nous pouvons donc : 

1. Obtenir une représentation explicite du polynôme d’interpolation : 

P n (x) = f[x 0 ] + 

n 

T. f [x 0 , Xi, • • • , x k ] (x - x 0 ) • • • (x - x k _i) 

k=l 

(Cette représentation est appelée représentation de Newton du polynôme) 

2. Calculer itérativement les valeurs P 0 (x),Pi(x), • • • ,P n (. x ) en un point x donné 


Si nous définissons 


f\Xi, Xi. |_i, ' • • , Xj 


S 

k=i 


fM 


(x k - Xi) • • • (x k - x k -i) (x k - x k+ i) ■ ■ ■ (x k - Xj) 


nous pouvons montrer que : 
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f ï,y, rp ry> 1 f [ X i + 1 i X i + 2 ' J X ! - :r / ■ I -' / ‘ I 7 I 1 

J [M , M • 1 • • • • ! •< iJ 

Cette relation suggère de calculer les différences divisées, à l’aide d’une table triangu- 
laire, de la façon suivante : 


Xi 

f(Xi) = f[Xi 

/[xi,a; i+ i] 

Xi+1, X^-i- 2 ] 

Xq 

f M 



X\ 

/N 

/[®o,a;i] 


X2 

/N 

/[a?i, aj 2 ] 

/[x 0 ,Xi,x 2 ] 


Exemple 


Déterminer le polynôme P 3 d’interpolation de Lagrange de / aux points ( Xi,f(xi )) 
suivants : (0,1), (1,3), (3, 2), (4, 5) 


Xi 

/M 




0 

1 




1 

3 

2 



3 

2 

1 

2 

5 

6 


4 

5 

3 

7 

6 

1 

2 


Le polynôme P 3 est donc donné par : 

P 3 (x) = 1 + 2 x — | x (x — 1) + 3 x (x — 1) (x — 3) 

Soit en développant l’expression de P 3 dans la base canonique : 

P 3 (x) = \x 2, — |ï 2 + yl + l 

Si on ajoute un point d’interpolation, il suffit de compléter le tableau des différences 
divisées 



P 4 (x) = P 3 (x) — | x (x — 1) (x — 3) (x — 4) 

— 5 ~,4 1 43 t 3 56 ~,2 i 43 ™ _ j_ i 

^ tv | 0 (Ay ^ tU | ^ tv | -L 

Corollaire 


- L’erreur d’interpolation de Lagrange est donnée par 
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e n (x) = /(x) - P n (x) = F(x) f[x 0 , ■■■ ,x n ,x] 

-Si /G C( n+1 )[a, b] alors il existe Ç x G [a, 6] tel que : 


f[x 0, • • • ,X n ,X 


/ (TC+ 1 ) fc) 

(n+1)! 


Preuve 


Pour x fixé, le polynôme Q(t) d’interpolation de Lagrange de / aux points x 0 , • • • , x n , x 
s’écrit : 


Q(t) = P n (t) + /[x 0 , ■■■ ,x n ,x](t-x 0 ) ■ ■ ■ (t - x n ) 


- Pour t = x, il vient 

f(x) = Q(x) = P n (x) + F(x) f[x 0 , ■■■ , x n , x] 

D’où la première relation 

- la deuxième relation du corollaire découle immédiatement de la première. 
Calcul de P n (a), a G [a, b] 

P n {oî) = f[x o] + f[x 0 , Xi] ( a - x 0 )+ 

f[xo, xi,x 2 \ (a - x 0 ) {a-xi) h 

f[xo,x i,-- - ,x n ] (a -x 0 ) (a -x^ (a -x n _i) 

Cette relation se réécrit de la manière suivante : 

P n (a) = f[x 0 \ + (a - x 0 ) (f[x 0 , xi] + {a - xi) (/[x 0 , xi, x 2 ] + 

■■■(a- x n _ 2 ) (/[x 0 , xi, • • • x n _i] + (a - x n _i) (/[x 0 , x x , • • • x n ]) 

•••))• 

L’évaluation se fait de droite à gauche de la manière suivante : 

ALGORITHME DE HÔRNER 

bo = f[xo,x i, • • -x n ] 
h = f[x 0 , x U ’” x n _i] + (a - x n _i) b 0 


bt f[x 0) Xi, • • • x n _j T {ol x n —i) bi—i 
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i — 1, • • • ,7i — l 


K = f[x o] + (a - x 0 ) b n - 1 = P n (a) 

Il est facile de vérifier que le coût opératoire pour calculer P n (a) est de : 

1. Algorithme de Newton : n(n+ 1) soustractions 

n( -" +1 ' > divisions 
n ( - n 2 +1 ' > divisions 

2. Algorithme de Hôrner : n additions 


n soustractions 


n multiplications 

Soit au total : n 2 + 3 n add/sous, n mult, div 
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CHAPITRE 3 


INTEGRATION NUMERIQUE 


Notations 

P„ : ensemble des polynômes à une variable, à coefficients réels, de degré inférieur ou 
égal à n. 

[a, b] : intervalle fini. 

/ : une fonction numérique définie et continue sur [a, b], 

{xo, x\, ■ ■ ■ , a; n } : n + 1 points distincts de [a, b], 
fl : une classe de fonctions définies sur [a, b]. 

Problème : Déterminer des formules de quadrature de la forme : 

/ b n 

f (x) dx = ^2 A i f( X i ) + E n(f ) 0 ) 

i=0 

n 

où Ai f(xi) est une valeur approchée de l’intégrale I. 

i= o 

E n (f ) l’erreur de quadrature associée 

avec les paramètres (Aj, Xi) 0<i<n calculés de manière que la formule (*) soit exacte sur 
fl, c’est-à-dire : 

V/ e fl , E n (f) = o 

3.1 Méthodes de Newton-Cotes 


1. Construction 
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Les noeuds {xi} 0<i<n étant choisis équidistants ( h = — , ay = a + i h , i = 
0, ••• , n). 

Remplaçons / par son polynôme d’interpolation de Lagrange p n (x) relativement 
aux points x 0 , xi, ■ ■ ■ , x n . 

De la relation : 


\/x G [a, b] , f(x ) = p n (x) + e n (x) 

où e„,(a;) représente l’erreur d’interpolation, il vient 

rb pb pb 

I — I f(x)dx— / p n (x) dx + / e n (x) dx (tHt) 

J a J a J a 

n 

On a : p n {x ) = f(xj ) Lj{x ) où {Li} 0<i<n désigne la base des polynômes de 

i = o 

Lagrange de P n . La relation se réécrit : 

/ b n 

f(x) dx = A i U) f( x i ) + E n(f) 
i=0 


2 . 


avec 


— I Li(x) dx , 


E n (f ) = / e n (x) dx 


0, 1, • • • , n et 


J a 

La formule (**) de quadrature est exacte sur P„ car si / G P n , alors f — p n 
en conséquence E n (f ) = 0. 

Calcul des poids 

Les coefficients ou poids A ^ de la formule de quadrature (**) sont données par le 

Théorème 


On prend : ay = a + i h , i — 0, 1, • • • , n avec h = — une subdivision de l’intervalle 
[a, b] en n sous- intervalles égaux. 


Vf = 0, 1, • • • , n on a : 


a (n) - A (n) 


i! (n— i)! 



j) du 
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Preuve. 


On a : A\'^ = / L^x) dx avec L^x) = JJ 3 

j: () j/‘ 

On en déduit que : 


OC 2 OC j 




n 


~n 1 / JJ ( X-Xj)dx 

JJ (Xi -Xj) “ j=0 jfr 


On a 


[Xi - Xj) = nij- xj 


) JJ {Xi-Xj) 


= (i! h 1 ) ((— l) n * ( n — i)\ h n ®) 
= (— l) n-i i\ (n — i)\ h n 


En faisant le changement de variable u = dans l’intégrale, on trouve : 

/ b n 

I l ( x — Xj) dx 

rn n 

= / JJ (a + u h ■— a — j h) h du 
° j=0 rfi 

/ n n 

| [ (u — j ) du 

•' "/// 

Pour établir la symétrie des coefficients, on fait le changement de variable : v — n—u 


pn n 

dans f JJ 

JO n 


( u — j ) du 


il vient 


rn n /*0 n 

/ n ( u ~j) du = - ji ( n — j — v ) dv 

J ~ 0 j¥=i 

çn n 

= (-!) n / JJ (■ v-j)dv 

J 0 — n 


0 
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d’où : = A ■ n ' 1 puisque (— l) n+î = (—1)" 1 


Exemple 1 : n — 1 Formule des Trapèzes 

A = A^ — h f u du = — 

./n ^ 


donc 

£4 (1| /( 


M = b -A [/(«)+ /w] 


i=0 


Exemple 2 : n = 2 Formule de Simpson 

4 2) = ^ 2 2) = | f {u-l)(u-2)du=^ 


A^ = —h I u{u — 2) du = — 

./n û 


4 h 


donc 

Ea <2, /( 


i=0 


Xi) = h —^~ [/(«) + 4 /(~^-“) + /(&)] 


3. Erreur d’intégration 


- Pour la formule des Trapèzes 

Pour évaluer le terme d’erreur pour la formule des Trapèzes, nous commençons 
par le lemme suivant 

lemme 

Soit / G C( 2 \[a, b]). Vx G [a, b] , 3£ æ G [a, b] tel que : 
f(x) = /(a) + (x - a) f{b) b ~J a {a) 

~ \{x-a ) ( b-x ) /” (Q 

De plus l’application x f”(Ç x ) es t continue sur [a, b]. 

Preuve. 

On pose : = f(x) — /(a) — (x — a) 

— c (x — a) (b — x) 
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Soit x 0 e]a, b[. 

On choisit c tel que : <L(a;o) = 0 

La fonction <f> G C^ 2 ^([a, b]) et s’annule en a,b,x o- 

Le théorème de Rolle montre que admet au moins 2 racines dans [a, b]. 

D’où la fonction <f>” admet au moins une racine dans [a, b]. Soit £ Xo cette racine. 
On a : 


0 = *”(&„) = /”(6j + 2c 

d’où C= -§/”(& o). 

Il est immédiat de voir que : x — > /” (£ x ) est continue pour tout x ^ a et x ^ b. 


Pour x — y a + 

i P fai - MzM 

lim - J F {F) = — y — /i»e 

x— >a+ 2 0 — O 

Pour a; — > 

1 -/'(&) + 

lim --f”fa) = ^L±± *=2- /ime 

i->t' 2 b — o 


D’après la règle de l’Hospital, on peut donc prolonger x — » F (F) P ar continuité 
en a et en 6. 


Formule des Trapèzes : 

si f G C (2 1([ a, b]) avec h = b — a 

y^f(x)dx = ^[f(a) + f(b)]-Pf' 2 >«), 

f e [a, b] 

Preuve. 

Ei (f) = J f( x ) dx - l ^~Y~ (f(a) + /(&)) 

= f (f( x ) - /(«) - (x - a) /(6) ~ /(o) ) d* 

= -.§ f f { 2 ) (F) i x - «) (b - x) do; 

J a 

La fonction x — » (x — a) (6 — x) a un signe constant sur [a, b] et x — » f^(F) es t 
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continue. En appliquant donc la formule de la moyenne sous forme intégrale, on 
trouve : 


E,(f) = — | / <2) (Ç) J (x-a)(b-x)dx 
= -àf m (0 (b -a? 


- Pour la formule de Simpson 

Un calcul plus compliqué que pour la formule des Trapèzes donne l’évaluation de 
l’erreur E 2 {f) 

£2(/) = -i/ (4) K) {b -a)* 

Formule de Simpson : 

sif eC^([a,b]) avec h = pp 

^ b f(x)dx = |[f(a) + 4 f(Up) + f(b)] - Tf(4) K)i 

t e [a, b] 

Preuve, (admise) 

7T 

Exemple : I = / sin(x) dx = 1 

J o 

Trapèzes : 

I = f (sm(0) + sin( f )) + sin(^) 

I = f + g sin(£) = 0.78539 + 0.32298 sin(£) 

e g [o, f ] 


Simpson : 

I = p (sm(0) + 4sm(|) + sm(|)) - sm(£) 
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I=u 1 + 2^2) - sâo «"(«) 

= 1.00227- 3.32052 x ÎO" 3 sin(Ç) 

e g [o, f ] 

3.2 Méthodes de Newton-Cotes composites 

Les formules de Newton-Cotes sont peu adaptées au calcul d’intégrales sur un grand 
intervalle d’intégration 

Exemple : 

e x dx = e 4 - e° = 53.598.. 

La valeur approchée obtenue par Simpson est : 

| (e° + 4 e 2 + e 4 ) = 56.769 ce qui est loin de la valeur exacte. 

On peut penser à augmenter n de telle façon à diminuer le pas d’intégration h = — 
mais cette solution s’avère numériquement mauvaise. 

En effet, pour les formules de Newton-Cotes de degré n : 

- Le calcul des devient difficile quand n augmente. 

- L’erreur de quadrature E n (f ) ne tend pas nécessairement vers 0, quand n tend vers 
l’infini pour une fonction donnée / G C([a,b\). 

11 est préférable de subdiviser l’intervalle d’intégration en un grand nombre de petits 
sous-intervalles, et d’appliquer sur chacun d’eux la méthode des Trapèzes ou de Simpson. 

On considère une subdivision de l’intervalle [a, b] en N intervalles égaux : 

U — a + ih, i — 0, 1, • • • , N et h — ^ 

On décompose : 

r-b N nti 

1= f(x)dx = f (x)dx 

J a i J ti — i 

pu 

Chaque intégrale / f{x)dx étant calculée par une formule de Newton-Cotes de bas 

J U - 1 



7 



Brahim Benouahmane 


degré. 


Formule des Trapèzes composites : 


si f G C^([a, b]) avec h = 
f f ( x )dx= ^[f(a) + 2 ^f(ti) + f(b)] 

•' a i=i 


b— a 
12 


h 2 f‘ 2 >«). 


{ e (a- b] 


Preuve. 

Pour i — 1, 2, • • • , N 

J f( X ) dx= \ (/(*<- O + /&)) - ^ / (2) (6)r 

Ci ^ [^i— lj^i] 


On en déduit 


h 


N - 1 


/(x) dx = - (/(a) + 2 ^ f(ti) + /(&)) 


2=1 


AT 


-§£/ <2) «.) 


2=1 


On voit facilement que : 

1 * 

min f (2 \x) < ^5I/ (2) (Ci) < 
æe[a,6] 'v ^ ' 


2=1 


max f {2) ( x ) 

x£[a,b\ 


Il existe donc £ G [a, 6] tel que : 

N 

/ (2> (?) = îE /^(£i) en appliquant à la fonction continue le théorème des 

i = i 

valeurs intermédiares. 

w L_ 

Il vient alors : g ^ / (2) (£i) = ^ / (2) (0 

2=1 

Remarque. 
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On appelle erreur la quantité : — ^ h 2 /^(£) 

| -^/r 2 / (2) (0 |< K h 2 

avec K = ^ (b — a) max | f^ 2 \x) \ 

x£[a,b\ 

On dit que la méthode des Trapèzes est d’ordre 2. 


Formule de Simpson composites : 


si f G C^([a, b]) avec h = et M = | 

„ b h M-l M 

/ f(x)dx= - [f (a) + 2 Y f(t 2 i)+4^f(t 2i _ 1 )+ 

•' a i=i i=i 

f ( b )]-Tg h4fl4| (î), îe[a,b] 


Preuve. 


Pour M = ^ avec N pair, on fait la décomposition 


"b M ' 

f(x) dx = / f(x) dx en M sous-intervalles de longueur 2 h. 

I i = \ Jt2i-2 


On en déduit 


h 


M 


f( x ) dx — — Yxm^) + 4 f(t 2 i-i) + f(t 2 i )) 
2=1 
M 

M-l M 

= f [/(°0 + 2 ^ f(fai) + 4^ /(Oî-l) + 


2=1 


2=1 


M 


m ]-!;£/ (4) fâ 


2=1 


La continuité de /-T entraîne l’existence de 
£ G [a, 6] tel que 


M 


i f E/ ,4) (« = / (4, (0 


2=1 


Le terme d’erreur s’écrit donc sous la forme 
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-SE / (4> «.) = ir / (4> «) = '* 4 / ,4, ( o 

Remarque. 

i-yfc 4 / < 4 ) K)i</i-fc 4 

avec A' = (6 — a) max j f^\x) \ 

x(z[a,b] 

On dit que la méthode de Simpson est d’ordre 4. 

Exemple : 

e x dx = e 4 — e° = 53.59815.. 

Déterminer l’entier N pour être sûr d’obtenir une erreur de quadrature inférieure à 
10 ~ 5 . 

Trapèzes : 

I T? h 2 / (2) (0 |= 3^ | | < J, e 4 < 10“ 5 

si N > 5397 

Simpson : 

I — h 4 1= 256 , I e ç I < 256 , e 4 < 10~ 5 

I 180 H J 60 I 45 V 4 I e I — 45 N 4 e — 1U 

si A^ > 76 

En fait, pour obtenir une précision absolue de 1CT 5 , il faut prendre N = 2674 pour la 
méthode des Trapèzes et A^ = 54 pour la méthode de Simpson. 
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CHAPITRE 4 


EQUATIONS NON-LINEAIRES 


Notations 

I : intervalle de M. 

[a, b] : intervalle fermé et borné de M. 

/ : une fonction numérique réelle définie et continue sur I. 

4.1 Séparation des racines 

Problème : Déterminer les racines de l’équation f(x ) = 0 sur / : 

(a e I , /(a) = 0) 

Une méthode de résolution numérique comportera deux étapes 

1. La séparation des racines ; cela consiste à déterminer un intervalle [a, b] C I tel que 
[a, b] contienne une racine a et une seule de l’équation f(x) = 0. 

2. Le calcul successif des racines séparées. 

La séparation des racines peut se faire essentiellement par : 

1. Une technique de balayage de / : 

- On découpe / en n intervalles [xi,x i+ 1 ] 

- On calcule successivement f(xi ) f(xi + \) ; si c’est strictement négatif, c’est qu’il 
existe un nombre impair des racines sur ]aq,aq + 1 [ 

- On recommence le découpage sur chaque intervalle où on a détecté au moins une 
racine jusqu’à "l’assurance" d’avoir isolé une racine. 
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2. ou graphiquement. 

Exemple 

f(x) = e x sin{x) — 1 I = [— n, n] 

L’intersection de fi(x) = sin(x ) et /^(x) = e~ x permet de déceler deux racines « 1 , 0:2 
avec 

«i e]0.5, 0.6[ et «2 e]3, 3.1[ 

4.2 Calcul d’une racine séparée 

On se place dans la situation suivante : 

Il existe un intervalle [a, b] C I tel que a soit la seule racine de f(x) = 0 sur ]a, b[ et 

/(«) /O) < 0. 


4.2.1 Méthode de Dichotomie (Algorithme) 

Entrée : £ (la précision désirée). 

Nq (le nombre maximum d’itérations) 

Sortie : valeur approchée de a ou un message d’échec. 

Poser n — 1 

Tant que n < N 0 et ^ > £ faire 

Poser p = 

Imprimer (n, a, b, p ) 

Poser n — n + 1 
Si /(a) f(p ) > 0 Alors 
Poser a = p 

Sinon 

Poser b = p 

Fsi 

Ftque 
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Si n > N 0 Alors 

Imprimer (la méthode a échouée après N 0 

itérations) 

Sinon 

Fsi 

Remarque 

Dans un algorithme on impose toujours un nombre maximal d’itérations afin d’éviter 
les boucles sans fin causées soit par une convergence trop lente soit par une erreur sur les 
données initiales. 


A chaque itération, l’algorithme construit un nouvel intervalle [a n ,b n \, autour de la 
racine cherchée a, qui est de longueur égale à la moitié de la longueur de l’intervalle 
précédent, c’est-à-dire : 

L ~ ^n — 1 ^n— 1 

Un U'n 2 

et donc on a : 


CLri. 


b— a 
2 n ~ 1 


Par conséquent, nous pouvons calculer le nombre maximal d’itérations pour obtenir 


bn 


2 — 


< £ 


On a : < £, d’où l’on tire : 

^ Log(b-a)~Log{e) 
11 ~ Log{ 2 ) 


Exemple 

La fonction f(x) = a; 3 + 4 x 2 — 10 a au moins une racine dans l’intervalle [1,2] puisque 
/( 1) = -5 et / (2) = 14. 


La valeur de a à 9 chiffres significatifs est 1.36523001. 
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Si on voulait obtenir ce résultat il faudrait faire 29 itérations. 

Il est intéressant de remarquer que l’approximation de a obtenue pour n — 9 est 
meilleure que pour n = 12. 


n 

Pn 

1 Pn - P29 | 

9 

1.36523438 

0.00000437 

12 

1.36499023 

0.00023978 

29 

1.36523001 

- 


Les avantages de cet algorithme sont la convergence assurée et une marge d’erreur sûre. 
Par contre il est relativement lent et il n’est pas toujours facile de localiser un intervalle 
sur lequel on a un changement de signe. Cet algorithme est souvent utilisé pour déterminer 
une approximation initiale à utiliser avec un algorithme plus rapide. 

4.2.2 Méthode de Newton-Raphson 

Soit x G [a, 6] une approximation de la racine a. On suppose donc que | x >- a | est 
petit et aussi que f'ix) 0. 

Si / G C^ 2 \[a,b]), le développement de Taylor-Lagrange de / à l’ordre 2 au point x 
s’écrit : 

f(x) = f(x) + (x-x) f'ix) + /”(Éx) 

où f x est un point entre x et x. 

En particulier, pour x — a, on obtient : 

/(«) = 0 = f(x) + (a - x) f'ix) + i -^r 1 /” (^ Q ) 

En négligeant l’inhniment petit (a — x) 2 d’ordre 2, on peut s’attendre à ce que : 


soit une meilleure approximation de a que x. 

Théorème 

On suppose que / G Ci 2 \[a,b]) avec a G]a, &[. 

Si oi est une racine simple de l’équation f{x) = 0 alors il existe un réel 9 > 0 tel que 
pour tout io G [a — 9, a + 9], l’itération de Newton : 
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x k+ i = x k ~ j&L k>0 
commençant en Xq converge vers a. 

Preuve. 

/'(a) 7^ 0 puisque a est une racine simple. 

Comme f est continue sur [a, b], il existe alors 6\ > 0 tel que : 

Va; G [a — a + Q\\ C [a, b] , f'(x) ^ 0 
Il suffit de montrer que la fonction g définie par 
g(x) =x-j^ r) 

définie et continue sur [a — 6\,a + ffi] vérifie les conditions suffisantes pour être une 
contraction. 

Il est facile de vérifier que g est dérivable sur [a — 0 1 , a + 9i], avec g' continue et définie 
par : 

n u \ _ fjx) r(x) 

Puisque g' {a) = 0, la continuité de g' assure l’existence d’un 6 > 0, 6 < 9\ tel que : 
Va; G [a — 6, a + 9\ , | g'[x ) |< K < 1 

Il reste à vérifier que g , définie sur [a — 9, a + 6\, est aussi à valeurs dans ce même 
intervalle 

Or Va; G [a — 9, a + 9] on a : 

| g{x) - a | = | g{ x) - g {a) | = | g'(Ç) | | x - a \ 

< K \ x — a \<\ x — ot\<6, £ est entre x et a. 

4.2.3 Algorithme de Newton-Raphson 

Entrée : Une approximation initiale p 

e (la précision désirée). 

N 0 (le nombre maximum d’itérations) 

Sortie ; valeur approchée de a ou un message d’échec. 
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Poser n — 1 

Tant que n < N 0 et | |> e faire 

Poser p = p - 
Imprimer ( n,p,f(p )) 

Poser n — n + 1 

Ftque 

Si n > N n Alors 

Imprimer (la méthode a échouée après N 0 

itérations) 

Sinon 

Fsi 


Exemple 


Considérons le même problème qu’en 2.1 
c’est-à-dire : f(x) = x 3 + Ax 2 — 10 


On a : f'(x) = 3 x 2 + 8 a; 


d’où : p n+ 1 = p n 


Pn+±Pn ~ l 0 


pour n = 5 le résultat est précis à 9 chiffres significatifs 


n 

Pn 

f{Pn ) 

1 

1.5000000000 

2.3750000000 

2 

1.3733333333 

0.1343454815 

3 

1.3652620149 

0.0005284612 

4 

1.3652300139 

0.0000000083 

5 

1.3652300134 

0.0000000000 
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4.2.4 Méthode de la sécante 

Il est parfois ennuyeux dans l’utilisation de la méthode de Newton- Raphson d’avoir à 
calculer f'(p„). L’algorithme suivant peut être considéré comme une approximation de la 
méthode de Newton-Raphson. 

Au lieu d’utiliser la tangente au point p n nous allons utiliser la sécante passant par les 
points d’abscisses p n et p n -i pour en déduire p n+ 1 . 

L’équation de la sécante s’écrit : 

s(x) = f(p n ) + (x - Pn) 

Si s(p n+ 1 ) = 0, on tire : 

Pn + 1 = Pn ~ f(pn)-7tPn- 1 ) 

4.2.5 Algorithme de la sécante 

Entrée : Deux approximations initiales p 0 et p\ 
e (la précision désirée). 

N 0 (le nombre maximum d’itérations) 

Sortie : valeur approchée de a ou un message d’échec. 

Poser n — 2 

<7o = f{po) 

Qi = f(Pi) 

Tant que n < N 0 + 1 et j q x \> e faire 

Poser p = pi — n ~ Po q, 
r n qi~qo 

Imprimer ( n,p,f(p )) 

Poser n — n + 1 
Po =Pi 

Pi = p 
Qo = Qi 
Qi = f(pi) 
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Ftque 

Si n > N 0 + 1 Alors 

Imprimer (la méthode a échouée après N 0 

itérations) 


Sinon 
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Remarque 

Si p n et p n -i sont proches (on peut espérer que ce sera le cas si la méthode converge), 
on peut considérer que 

l) 

Pn Pn—1 

est une approximation de f{p n )- Dans ce cas l’algorithme de la séante est une ap- 
proximation de l’algorithme de Newton. 

Exemple 

Considérons le même problème que précédemment, 

c’est-à-dire : f(x) = x 3 + Ax 2 — 10, avec comme point de départ p 0 = 2 et pi = 1 

Nous constaterons que la convergence est moins bonne que celle obtenue par la mé- 
thode de Newton. 


4.2.6 Méthode du point fixe 

On remplace la résolution de l’équation 
f(x) = 0 (E,) 

par la résolution de l’équation 
g(x) = x (. E 2 ) 

Le choix de g doit être tel que : 

1. la suite itérative définie par : 

x 0 G [a, b] 

x n+l = g(x n ) n > 0 
soit convergente. 

Si g est continue la limite a de la suite (x n ) n >o est alors solution de (E 2 ). 

2. la solution de (E 2 ) soit aussi solution de ( E 1 ). 

Nous pouvons observer que la méthode de Newton s’écrit toujours : 
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Pn + 1 = g(Pn ) OÙ 5f(x) = X ~ 

Si la fonction g est continue et si l’algorithme converge, c’est-à-dire si p n tend vers a, 
alors a vérifie : a = g (a). 

On dit que a est un point fixe. 

Théorème 

On considère une fonction g définie et continue sur [a, b], à valeurs dans [a, b]. On 
suppose que g est une contraction sur [a, b] : c’est-à-dire qu’il existe K G [0, 1 [ tel que 

V(x, x ') G [a, b] x [a, b] : 

| g(x) — g(x') | < K \ x — x' \ 

Alors l’équation (E 2 ) admet une racine unique a G [a, b] et on a 
\/x 0 G [a, b] la suite ( x n ) n > 0 définie par 
X-n+l Çl^Xn) 
converge vers a et 

\a-x n \< \ x 1 -x 0 \ 

Remarque 

Ce théorème est un résultat de convergence globale car la convergence de la suite 
{x n ) n > o vers a est assurée quel que soit l’initialisation x 0 . 

Preuve. 

Montrons que la suite définie dans i) converge vers une racine de (E 2 ). 

Vj > 0 | x j+1 - Xj | = | g{xj) - g(xj- 1 ) | 

< K | Xj — Xj _ i |< JO | xi — x 0 | 
pour p > 1 fixé, il vient : 

n+p— 1 

| %n+p %n | ^ ^ ^ | Œj+l Œj | 

j—n 

p— 1 oo 

<| Xi - x 0 I K" Rj <| an — a?o | K" Y Rj 

3=0 3=0 
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soit | x n+p -x n \< ^ | Xi - x 0 | (*) 

On en déduit que : lim | x n+p — x n |= 0, c’est-à-dire que (a: n ) n >o est une suite de 

n— ï-+oo 

cauchy dans M et par suite convergente vers a G [a, b\. 

De la continuité de g, il découle que : 
a = lim x n+1 = lim g(x n ) 

n— >•+ oo n — >-+oo 

= g( lim x n ) = g (a) 

n— >+oo 

Montrons maintenant que (E 2 ) possède une racine sur [a, b], suposons l’existence de 
deux racines distinctes a et a', alors 

| a — a' |=| g (a) — g(a!) \< K \ a — a' \ 

ce qui aboutit à la contradiction 

(1 — K ) j a — a! |< 0 

En faisant tendre p vers 00 dans l’inégalité (*) on trouve 
\a-x n \< \xi — x 0 \ 

Corollaire 

On considère une fonction g définie et continue sur [a, b], à valeurs dans [a, b]. 
g est une contraction si les conditions suivantes sont vérifiées : 

1. g est dérivable sur |a,b] 

2. Va: G [a, b] , | g'{x ) |< K < 1 

Preuve 

De la formule des accroissements finis, il vient que 
\/(x,x r ) G [a, b] x [a, b], il existe £ Gjaqa/f tel que : 
g(x) — g(x') = g'(0 (x - x') 
d’où : \/(x,x') G [a, b] x [a, b], 

| g(x) — g(x') |=| #'(£) \ \ x — x' \< K \ x — x' |. 
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4.2.7 Algorithme du point fixe 

But : Trouver une solution de x = g(x) 

Entrée : Une approximation initiale p 

e (la précision désirée). 

Nq (le nombre maximum d’itérations) 
Sortie : valeur approchée de a ou un message d’échec. 

Poser n — 1 

Tant que n < N 0 et | g(p) — p \> £ faire 
Poser p — g(p) 

Imprimer (n, p , g(pj) 

Poser n — n + 1 

Ftque 

Si n > Nn Alors 

Imprimer (la méthode a échouée après N 0 

itérations) 

Sinon 

Fsi 


4.2.8 Exemple 1 

f(x) = x 3 + 4 x 2 — 10 , [a, b] = [1, 2] 

On considère : 

gi(x) = 10 + x — 4 x 2 — x 3 
g 2 (x) = -Ax 
g 3 (x) = | \/10 - x 3 

9±{x) = ^ 
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On pourra vérifier que g 3 et g 4 sont des bons choix pour résoudre l’équation f(x) = 0 
ce qui n’est pas le cas de g\ et r/ 2 . 

Exemple2 

f(x) = x 2 — 2x — 3 

L’équation f(x ) = 0 admet deux racines qui sont -1 et 3. 

On considère : 

g(x) = \/2x + 3 

On utilise un estimé Xo = 4, on obtient : 

X! = Vïï ~ 3.316 
a ; 2 = V9.632 ~ 3.104 
£3 = x/9.208 ~ 3.034 
x 4 = 79.068 ~ 3.011 
x 5 = 79.022 ~ 3.004 

(x n )n> 0 apparaît comme une suite convergeant vers 3. 

On considère : 

ff( x ) = ^2 

Si x 0 = 4, alors : 

X\ = 1.5 
x 2 = -6 
0:3 = -0.375 
x 4 = —1.263 
X 5 = —0.919 
x 6 = —1.028 
x 7 = —0.991 
x 5 = —1.003 
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On voit que la suite ( x n ) n > 0 converge vers -1. 

On considère : 

g(x) = ^ 

Pour x 0 = 4, on obtient : 

X! = 6.5 
rr 2 = 19.635 
x 3 = 191.0 

Cette suite diverge. 

2.2 Remarque 

Il est parfois difficile de transformer f(x) = 0 en x = g(x) par des transformations 
algébriques simples. 

Dans ces cas on pourra réécrire f(x) = 0 comme 
x = x + k f(x ) = g(x) 

où k est une constante non nulle que l’on choisit 

On choisit k afin que la dérivée de g(x) soit telle que : 

| g\x) | = | 1 + kf'(x) \< 1 

afin d’assurer la convergence. 

(Prendre k = \ pour l’exemple2) 
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CHAPITRE 5 


RESOLUTION NUMERIQUE DE PROBLEMES DIFFERENTIELS 


5.1 Problème de Cauchy 

Le but de ce chapitre est d’étudier des méthodes pour approximer les solutions y(x) 
d’un problème de type 

£ = f(x,y) (i) 

y(x 0 ) = y 0 (2) 

où la fonction f(x,y) est une fonction connue de deux variables. 

(1) détermine une famille de solutions y = <h(a;, c) (appelées combes intégrales) dé- 
pendant d’un paramètre. 

(2) permet de choisir le membre de cette famille de courbe passant par (x'o, yo) (ceci, 
en déterminant une valeur particulière de c). 

Du point de vue géométrique, l’équation (1) signifie que la pente de la courbe intégrale 
y = <h(a:,c) passant par un point (x,y) quelconque est donnée par f(x,y). 


Définition 

On dit que la fonction / : [a, b] x I — » M est lipschitzienne en y dans [a, 6] x M s’il 
existe A > 0 tel que : 

I f(x,y) - f(x,z) |< A | y - z | 

Vx G [a, b] Vy, z G I 

A est appelée constante de Lipschitz. 
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Théorème (admis) 

Soit le problème de Cauchy ((1),(2)) 

Si / : [a, b] x M — y M vérifie les hypothèses : 

1. f continue 

2. f lipschitzienne en y 

alors le problème de Cauchy ((1),(2)) admet une solution unique. 

5.2 Approche numérique 

On prend : ay = a + i h , i — 0, > • • , N avec h = une subdivision de l’intervalle 
[a, b] en N intervalles égaux. 

On cherche N nombres 2/i , 2 / 2 , • • • ,Vn où y,; est une valeur approchée de y(xi). 

Puis on reliera ces points par interpolation pour définir une fonction y h sur [a, b]. 

On va voir dans ce qui suit trois méthodes numériques permettant de calculer la 
solution approchée y l+ 1 au point x i+ i en utilisant celle au point ay. 

5.2.1 Méthode d’Euler 

On considère que sur le petit intervalle 
[a; 0 , x 0 + h] la courbe n’est pas très éloignée de sa tangente en x 0 

z(x) = y(x 0 ) + y'{x 0 ) (: x - x 0 ) 

Une bonne approximation de y(x 1 ) est donnée par 
yi = Vo + hf(x 0 ,y 0 ) 

On considère ensuite que f(x\,yi) est une approximation de y'{x 1 ) et sur l’intervalle 
[x\,X2 ], on remplace la courbe par sa tangente approchée en X\ 

z(x) = y(x 1 ) + y'{x 1 ) (x - Xi) 

-yi + /(®i,î/i) (x-X!) 
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Une bonne approximation de y(x 2 ) est donnée par 
V2 = Vi + hf(x 1 ,y 1 ) 

A la ième étape en approximant y(xi) — yi, on obtient : 
y(xi + h) ~ y.i + hf(xi,yi) 
ce qui suggère de nouveau l’équation : 

Vi+i = Vi + hf(xi,yi) (3) 

appelée équation aux différences pour la méthode 
d’Euler. 


Exemple 


y' = —y + x + 1 

y(x 0 ) = 1 

La solution exacte est donnée par : 
y(x) = x + e~ x 

Prenons n — 10 , h = 0.1, (3) devient : 


Ui + 1 = 0.9 yi + 0.1 Xi + 0.1 i = 0, 1, • • • ,9 


On obtient alors le tableau : 


i 

Xi 

y(xi) 

Ui 

1 y{xi) - yi \ 

0 

0 

1.004837 

1.000000 

0.004837 

1 

0.1 

1.018731 

1.010000 

0.008731 

2 

0.2 

1.040818 

1.029000 

0.011818 

3 

0.3 

1.070302 

1.056100 

0.014220 

4 

0.4 

1.106531 

1.090490 

0.016041 

5 

0.5 

1.148812 

1.131441 

0.017371 

6 

0.6 

1.196585 

1.178297 

0.018288 

7 

0.7 

1.249329 

1.230467 

0.018862 

8 

0.8 

1.306570 

1.287420 

0.019149 

9 

0.9 

1.367879 

1.348678 

0.019201 


On peut remarquer que l’erreur croît légèrement lorsque x t croît. 
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Convergence 

Théorème 

Si / vérifie les hypothèses : 

1. / G C\[a,b] x R) 

2. f est lipschitzienne en y : 

I f(x,y) -f(x,z) |< K | y -z | 

Vx G [a, b] Vy, z G R , K > 0 

alors la méthode d’Euler converge. 

Plus précisément, si on pose : 

M — M ax{ | y'(t) |, t G [a, 6]}, 
on a la majoration : 

I e n | = | y n - y(x n ) \< eK(b ^ ] - 1 M h 
et lim Max{ I e n |, n = 1, • • • , = 0 

>-+oo 
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Remarques 

1. Le résultat du théorème précédent s’exprime sous la forme : 

| e n | < Ah , A > 0 une constante 
c’est-à-dire que la méthode d’Euler est du premier ordre. 

2. Une méthode d’orde 1 ne converge pas assez vite pour donner des résultats pratiques 
intéressants. 

Exemple 

y' = ly + x 2 e x , y{ 1) = 0 , x G [1, 2] 

On veut déterminer le pas h pour que l’erreur soit inférieure à 0.1. 

Puisque | f(x, y) - f(x , z) | = | l(y- z) \ 

Mxe [1,2], K = 2. 

D’autre part en résolvant l’équation on obtient : 
y(x) = x 2 e x — e x 2 , 
d’où : M = y\ 2) = 4 e (—1 + 2 e) 

La borne devient donc : 

e^_ 4 e (—1+2 e) , 

2 2 a 

Il suffira donc que h < g3 — 0.0011222 

5.2.2 Méthodes de Taylor 

Line première façon d’améliorer la méthode d’Euler (au sens où l’erreur variera au 
moins comme h 2 ) consiste à utiliser un développement de Taylor jusqu’à l’ordre 2. 

Si y(xi) est donné, on a : 

y{x i+ 1 ) = y(xi) + (x i+ i - Xi ) y'(xi)+ 

l (x i+1 - Xi) 2 y {2 \xi) + | (x i+1 - Xi) 2 y {3) ^i) 
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où & G [xi,x i+ 1] 

Si h = x i+ i — Xi est assez petit, on a donc l’égalité approchée 

y(x i+ 1) ^ y(xi) + hy'ixi) + Çz/ (2) (^) 

Puisque y'{x ) = f(x,y(x)) on a : 

y”(®) = ||(®,î/(®)) + ^(x,y(x))y'(x) 

= ^(x,y(x)) + | ${x,y(x))f(x,y(x)) 

Tout ceci suggère la méthode aux différences suivantes 
Yi+i =yi + hf(x i5 yi) + 

Ç [ff(xi,yi) + ||( x i,yi) f ( x i,yi)] (4) 
appelée méthode de Taylor d’ordre 2 
Remarques 

1. On peut montrer que pour chaque Xi fixé, l’erreur varie comme h 2 . 

2. Cette méthode possède le désavantage qu’elle nécessite le calcul de 

3. Cette méthode se généralise facilement, cependant les méthodes de Taylor nécessi- 
teront le calcul de 44 44 44 , £ / , 1 - 4 ., • • • , ce qui peut être très laborieux. 

ox 1 oy^ox zl oxoy 1 o y z 7 7 ^ ^ 


Exemple 

y 7 = -y + x + 1 
yo = y( x o) = i 

Puisque ^ = 1, = — 1 l’équation aux différences s’écrira : 

Vi + 1 = IJi + h (- y.i + x* + 1) + ^ [1 r- (-y* + x* + 1)] 
ï/i+l = Vi + h [(| — 1) 2/i + (1 — |) Xj + 1] 
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Prenons h — 0.1 on obteint alors le tableau suivant : 


Ul 

erreur 

1.0190251 

1.07082 

1.149404 

1.2949976 

2.942 x 10 -4 
4.820 x 10~ 4 
5.924 x 10~ 4 
6.470 x 10" 4 


5.2.3 Méthode de Runge-Kutta 

On voudrait conserver l’avantage des méthodes d’ordre supérieur mais corriger les 
inconvénients dûs au calcul des dérivées partielles de /. 

Ces différentes méthodes sont basées sur la formule de Taylor à plusieurs variables : 

f(x, y) = f(x o, y 0 ) + [(z -xq )% i + (y- yo) ]+ 
è [0 - x o) 2 fÿ + 2 (x - x o) (y - y 0 ) 

n + 1 

( X - X o) 2 0] + • ' ' + 0ï)T Y, C n + 1 K X - X o) n+1 ~ j 

3=0 

(: V - VoY K, ri)] (5) 

où ( € (x,xo) , ri € (y, yo) 

Parmi toutes ces méthodes, la plus utilisée est celle d’ordre 4 mais les calculs menant 
à l’algorithme sont très laborieux. 

Illustrons l’idée pour la méthode de Runge-Kutta d’ordre 2. 

L’idée est d’essayer de remplacer le terme 

f(xi,Vi) + \ [^i(xi,Vi) + ^(xi,yi) f(xi,yi)] 
dans la formule (4) par une expression de type 
af(xi + a,yi + p) 
en utilisant (5) on obtient : 
a f (xi + a,yi + (5) = af(x h yi) + a a fyO*, yi)+ 

a P | fai,yt) +aR 2 (Ç i ,ri i ) 
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avec & G (a*, x t + a) , ^ G (j/*, y t + /3) 

En identifiant les coefficients de / , ^ ^ on voit que l’on doit choisir : 

a = l , aa=\ , a P = § f(xi,Ui) 
c’est-à-dire : 

a = l , a = | , /3 = f f(xi,yi) 

En substituant dans (4) on obtient donc la formule aux différences : 
yi +1 = yi + hf(x; + |,yi + |f(xi,yi)) 


d’où l’algorithme suivant : 
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Méthode de Runge-Kutta d’ordre 2 

Soit N = XN ~ X ° =nornbre de pas nécessaire 
Poser i — 0 

Tant que i < N faire 

k = h f (xi, Di) 

Vi+i = Vi + h f(x.i + §, y t + |) 

Xi + 1 = Xi + h 

i = i + 1 

Ftque 

Ceci nous amène à la méthode de Runge-Kutta d’ordre 4 pour laquelle nous donnons 
directement l’algorithme : 

Méthode de Runge-Kutta d’ordre 4 

Soit N = XN ~ X ° =nombre de pas nécessaire 
Poser i = 0 

Tant que i < N faire 

h = h f(xi, yi) 
k 2 = hf(xi + | ,Vi + l h) 
h = hf(xi + | ,Vi + l k 2 ) 
h = h f(xi + h, y.i + k 3 ) 

Vi+i — Ui + | -f 2 k 2 + 2 k 3 + k^) 
x i+ 1 = Xi + h 
i = i + 1 

Ftque 

Commentaire 

Les méthodes de Runge-Kutta sont faciles à utiliser, cependant elles sont relativement 
coûteuses car elles nécessiteront plusieurs évaluations de f(x,y ) à chaque pas. 

Exemple 
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1. y' = — y + x 2 + 1 , 0 < x < 1 , 2 /( 0 ) = 1 


Solution exacte est : y(x) = — 2 e x + x 2 — 2 a; + 3. 


Utiliser la méthode de Runge-Kutta d’ordre 2 avec h = 0.1 , N — = 10 


Xi 

y(xi) 

Ui 

erreur : y* — y l 

0 

1.00000 

1.00000 

0 

0.1 

1.00033 

1.00025 

0.0000751639 

0.2 

1.00254 

1.00243 

0.0001122438 

0.3 

1.00836 

1.00825 

0.0001178024 

0.4 

1.01936 

1.01926 

0.0000974985 

0.5 

1.03694 

1.03688 

0.0000562001 

0.6 

1.06238 

1.06238 

0.0000019171 

0.7 

1.09683 

1.09690 

0.0000732812 

0.8 

1.14134 

1.14150 

0.0001548478 

0.9 

1.19686 

1.19710 

0.0002440317 

1.0 

1.26424 

1.26458 

0.0003386469 


2. y' = — y + x + l , 0 < rr < 1 , y( 0) = 1 


Solution exacte est : y(x) = —2 e~ x + x 2 — 2 x + 3. 

Utiliser la méthode de Runge-Kutta d’ordre 4 avec h — 0.1 , N = Ry = 10 

On peut comparer les résultats avec ceux obtenus précédemment pour voir l’aug- 
mentation de la précision. 


Xi 

y Ou) 

y* 

erreur : y* — y* 

0 

1.0000000000 

1.0000000000 

0.0000000000 

0.1 

1.0048374180 

1.0048375000 

0.0000000820 

0.2 

1.0187307531 

1.0187309014 

0.0000001483 

0.3 

1.0408182207 

1.0408184220 

0.0000002013 

0.4 

1.0703200460 

1.0703202889 

0.0000002429 

0.5 

1.1065306597 

1.1065309344 

0.0000002747 

0.6 

1.1488116361 

1.1488119344 

0.0000002983 

0.7 

1.1965853038 

1.1965856187 

0.0000003149 

0.8 

1.2493289641 

1.2493292897 

0.0000003256 

0.9 

1.3065696597 

1.3065699912 

0.0000003315 

1.0 

1.3678794412 

1.3678797744 

0.0000003332 
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RESOLUTION DES SYSTEMES LINEAIRES 


Introduction 


La résolution de grands systèmes (linéaires ou non-linéaires) est pratique courante de 
nos jours, spécialement en génie mécanique, en génie électrique, et de façon générale, dans 
tous les domaines où l’on s’intéresse à la résolution numérique d’équations aux dérivées 
partielles. 

Rappel sur les systèmes linéaires 

Un système de n équations à n inconnues peut toujours s’écrire sous la forme : 

Ax = b 

où A = (cii t j)i<i tJ < n est une matrice n x n 
et x et b sont des vecteurs colonnes de dimension n. 

Méthode de résolution 

La méthode de résolution la plus étudiée (et une des plus utilisées) s’appelle méthode d’élimination 
de Gauss . Elle a pour but de remplacer le système Ax = b par un système triangulaire 
supérieur avec des 1 dans la diagonale. C’est-à-dire de la forme : 


( 1 01,2 • • • Ôi , n ^ 


( Xi \ 


( h \ 

1 Ô 2 ,3 • • 0>2,n 


x 2 


î>2 

1 Cl n —i tn 


•^n— 1 


bn— 1 

V 1 ) 


\ X n ) 


{ K J 


On en déduit alors facilement les composantes du vecteur solution x : 
■t' n b n 

n 

Xi = bi — ^2 âi,j x j pour i — n — 1, • • • , 1 
j=i + 1 

Elimination de Gauss sur un exemple 


1 



! xi + x 2 = 3 
2 xi + x 2 - x 3 = 1 
3 xi - x 2 - x 3 = -2 

qui s’écrit : 



• Elimination de la première inconnue 


Eliminons X\ dans les deux équations 2 et 3. Pour ce faire on multiplie la première 
équation par et on soustrait à ces deux équations 

11 0 W Xl \ ( 3 

0 -1 -1 \ \ x 2 = -5 

0-4-1 ) \x 3 ) V - 11 

• Elimination de la deuxième inconnue 

On divise la deuxième équation par le 2 eme pivot 

11 0 W Xl \ f 3 

0 1 1 \ \ x 2 = 5 

0 -4 -1 J \ x 3 ) \ -11 

Eliminons X 2 dans la dernière équation. Puis divisons l’équation ainsi obtenue par 
le 3 ème pivot 


1 1 0 W xi \ / 3 

0 1 1 \ \ x 2 = 5 

0 0 1 ) \ x 3 ) V 3 

Finalement on trouve : 


X 3 = 3 x 2 = 2 x\ = 1 
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Algorithme : Elimination de Gauss 


Pour i — 1 jusqu’à n — 1 faire 
Pour j — i + 1 jusqu’à n faire 

,7 ~ „ . 

Fpour 
bi <- 

Pour k = i + 1 jusqu’à n faire 
Pour j =2+1 jusqu’à n faire 


^k,j ' Q'kJ @‘k,i ^i,j 

Fpour 

bk t— bk — a>k,i bi 

Fpour 


Fpour 


Cette première phase de la méthode de Gauss est la triangularisation. Il reste main- 
tenant à résoudre le système triangulaire supérieur obtenu après la triangularisation. 

x n t b n 

Pour i = n — 1 jusqu’à 1 pas -1 faire 

n 

Xi i bi ^ ^ C/.y Xj 
j=i + 1 

Fpour 

Notons toutefois que notre algorithme ne peut être exécuté jusqu’à la fin que si les 
pivots successifs sont tous non nuis. 

Nous pouvons nous demander s’il existe une relation entre la matrice de départ et les 
matrices triangulaires obtenues. Ce lien existe : 

Si U et L désignent les matrices triangulaires obtenues après triangularisation : 
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^1,1 

h,i 

h, 2 


n— 1,1 

ln- 1,2 • • 

ln— l,n— 1 

ln, 1 

ln, 2 

ln,n— 1 1 

/ 1 

«1,2 

• • ^l,n 


1 «2,3 

• • ^2 ,n 



1 Un— l,n 


V 1 J 

(et si l’algorithme d’élimination n’exige pas d’échange de lignes), on a : 


A = LU 

On dit dans ce cas que l’ona décomposé (ou factorisé) A en un produit L U. 

Nous ne démontrerons pas cette proposition. Nous nous contenterons de la vérifier sur 
notre exemple : 


( 

i 

0 

0 

0 



( 1 

1 

0 

3 



( 

1 

1 

0 

3 \ 


2 

-1 

0 

0 



0 

1 

1 

5 




2 

1 

-1 

1 


3 

-4 

3 

0 



0 

0 

1 

13 

3 




3 

-1 -1 

2 

V 

-1 

3 

0 

-13 

) 


1 0 

0 

0 

1 

) 


\ 

-1 

2 

3 

-1 / 


Il y a une classe importante de matrices pour lesquelles l’élimination peut toujours 
s’opérer sans échange de lignes (c’est-à-dire le pivot a hJ n’est jamais nul pendant l’algorithme 
d’élimination). Ce sont les matrices à diagonale strictement dominante. 

Définition 

Une matrice A est dite à diagonale strictement dominante si pour i = 1,2 
l’inégalité stricte 


est vérifiée. 


| ^ 'y ) «i y 
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Nombre d’opération pour l’élimination 
de Gauss 

Pour évaluer la rapidité d’un algorithme, il est important de connaître le nombre 
d’opération qu’il nécessite. 

On montre que l’algorithme d’élimination de Gauss fait appel à : 

• " ( ~ n 3 ~ 1 ^ additions et autant de multiplications 

• divisions 

La résolution du système triangulaire supérieur, quant à elle, nécessite additions 

et autant de multiplications. 

Au total la résolution d’un système de n équations linéaires par la méthode de Gauss 
demande : 

• " ( -' i '~ 1 g 2rt + 5 ) additions et autant de multiplications 

• divisions 

Remarques 

1. On voit que, pour n grand, le nombre d’additions et de multiplications est de l’ordre 

3 

de î |-. Ainsi si l’on multiplie par 2 la dimension d’un système il faudra 8 fois plus 
de temps pour le résoudre. 

2. Si l’élément diagonal a*,* est nul, on cherche, dans la même colonne, un élément 
d’indice plus grand non nul, puis on échange les lignes correspondantes. Si ceci est 
impossible, le système est singulier. 

3. On est parfois amené, pour des raisons de stabilité numérique, à effectuer des 
échanges de lignes même si a iti est non nul. Ceci conduit à des stratégies dites 
de pivots. 

Elimination de Gauss avec changement 
de pivots 
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Exemple 


On veut résoudre : 

0 1 3\/n\ ( l 

5 2 3 x 2 = 4 

6 8 1 ) \ x 3 ) V 1 

Echange de lignes 

Dans ce cas on ne peut pas directement appliquer l’algorithme d’élimination de Gauss 
puisque le 1 er pivot est nul. 

Toutefois, on peut échanger deux lignes pour obtenir un premier pivot non nul. 

On choisit comme premier pivot le plus grand coefficient de la l ere colonne. Pour cela 
on échange la l ere et la 3 eme ligne. 

6 8 1 \ / Xi \ fl 

5 2 3 x 2 = 4 

0 1 3 ) \ x 3 ) V 1 

Maintenant on peut appliquer l’algorithme d’élimination de Gauss 

• Elimination de la première inconnue 

La première opération consiste à diviser cette équation par a^i = 6 encore appelé 
premier pivot 


/ 1 f i 

L 3 6 

5 2 3 
V 0 1 3 


X\ 

X 2 

%3 



Eliminons l’inconnue X\ dans les deux équations 2 et 3. Pour ce faire on multiplie 
la première équation par a, l: \ et on soustrait à ces deux équations 
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• Elimination de la deuxième inconnue 


On divise la deuxième équation par le 2 eme pivot 



Eliminons x 2 dans la dernière équation. Puis divisons l’équation ainsi obtenue par 
le 3 ème pivot 


Finalement on trouve : 


1 


4 
3 

0 1 
0 0 



1 

6 

_ 19 

28 

47 

97 


47 44 67 

Xs ~ 97 X2 ~ ~97 Xl ~ 97 


Algorithme : Elimination de Gauss avec 
choix du pivot 

Pour i — 1 jusqu’à n — 1 faire 

Pivotage partiel 


Pour j — i + l jusqu’à n faire 


® i,j t 


a i,j 

0 * 1,1 


Fpour 



Pour k — i + 1 jusqu’à n faire 
Pour j = i + 1 jusqu’à n faire 

Fpour 
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Fpour 


Fpour 
b n < — — 

Ûn,n 

La procédure du pivotage partiel, introduite dans la l ere ligne de l’algorithme précédent 
s’écrit : 

k <— i 

m <r- abs(di t i) 

Pour j — i + 1 jusqu’à n faire 
s abs(aj y i) 
sf m < s alors 

k 4— j et m 4— s 
sinon 
Fsi 
Fpour 

Si k ^ i alors 

Pour j = i jusqu’à n faire 

t i — Ojj 

j ' 

dk.j t 
Fpour 
t <— bj 
bi <— bk 
bk ~ t 
sinon 


Fsi 



Factorisation 


Dans le cas sans pivotage, si l’on doit résoudre un système où le membre de droite 
change, il y a intérêt à effectuer la réduction à la forme triangulaire une fois pour toutes. 

En effet, si A = LU on peut résoudre : A x — b en résolvant : 

L z = b 
U x = z 

Les systèmes étant triangulaires, la résolution ne nécessite que l’exécution d’une re- 
montée et d’une descente triangulaire. 

Remarque 

Lorsque la matrice A est symétrique définie positive il y a intérêt à utiliser la méthode de Cholesky 
qui est une variante de la méthode de Gauss. 

La méthode de Cholesky consiste à effectuer 
la décomposition A = L L*, c’est-à-dire que 

U = LL On effectue alors deux fois moins d’opérations arithmétiques qu’avec la méthode 
de Gauss. 

Méthode de Crout 

Nous allons montrer sur un cas particulier que, l’on peut parfois factoriser directement 
une matrice A, en tenant compte de certaines structures particulières. 

Soit 

« 1,1 oq, 2 

^ 2,1 « 2,2 « 2,3 

«3,2 «3,3 «3,4 

«n— 1,77— 2 «71—1,71—1 «7i—l,7i 

«71,77—1 «71,77 

une matrice tridiagonale. Nous cherchons des facteurs L et U de la forme : 
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( kl 

ki 


h, 2 
h, 2 


\ 


L = 


13,3 


ln—l,n—2 ^n— l,7i— 1 

k,n— 1 k,n ) 


U = 


( 1 « 1,2 

1 «2,3 

1 «3,4 


\ 


1 «n— l,n 

1 / 


La multiplication matricielle conduit aux équations: 


1. /pi = ai,! (ligne 1 par colonne 1) 

2. Ijj-i Uj - ij + Ijj = djj (ligne j par colonne j) 

3. lj+ij = cij+ij (ligne j+1 par colonne j) 

4- lj,j «jj+i = a j,j + 1 (ligne j par colonne j+1) 

Ceci suggère l’approche suivante : 

Pour j — 1, 2, • • • , n 

• Déterminer la colonne j de L à l’aide de (b) 

((a) si j=l) puis de (c). 

• Déterminer la rangée j de U k l’aide de (d). 

Si on se réfère à (c) on voit que la sous-diagonale de L coincide avec celle de A. Il 
suffit donc de stocker la diagonale de L et la sur-diagonale de U, ce que l’on fait, bien 
sûr, dans deux vecteurs. 
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